Floppiness, cutting, and freezing: Dynamic critical scaling near isostaticity 
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The isostatic state plays a central role in organizing the response of many amorphous materials. 
We demonstrate the existence of a dynamic critical length scale in nearly isostatic spring networks 
that is valid both above and below isostaticity and at finite frequencies, and use scaling arguments 
to relate the length scale to viscoelastic response. We predict theoretically and verify numerically 
how proximity to isostaticity controls the viscosity, shear modulus, and creep of random networks. 
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Random networks of springs display two distinct 
phases distinguished by the presence or absence of floppy 
modes - zero frequency motions that neither compress 
nor stretch springs. Floppiness is related to network 
structure via the mean coordination z [T] ; networks with 
floppy modes reside below the isostatic coordination z c . 
Many amorphous materials possess an isostatic state, in- 
cluding fiber networks, covalent glasses, foams, and emul- 
sions. Hence the viscoelasticity of damped random net- 
works, which remains poorly understood, could poten- 
tially provide insight into a broad class of materials, in- 
cluding how structure relates to response [TUT^]. 

The dramatic impact of floppiness on response is ap- 
parent in numerical measurements of the creep compli- 
ance J(s) = 7(s)/ito, shown in Fig. 1; j(s) is the Laplace- 
transformed shear strain accrued after a small step stress 
<7o- At long times or vanishing s, hyperstatic networks 
(z > z c ) approach constant strain, J ~ s , while hypo- 
static networks (z < z c ), approach constant strain rate, 
J ~ 1/s. Thus networks with floppy modes are fluids, 
and those without arc solids. Moreover, networks' elas- 
ticity and viscosity clearly vary with proximity to z c . 

Together with prior work [1 \5\ HI ED EE2] , data such 
as in Fig. 1 strongly suggest that the isostatic state is a 
nonequilibrium critical point. We will show this is indeed 
the case, and that proximity to z c controls viscoelasticity 
in random spring networks. To fully understand dynam- 
ics near a critical point, it is essential to identify the 
associated diverging length scale [13] • While hyperstatic 
networks possess an "isostatic length" £* ~ 1/Az, where 
Az = z — z c 14-18], its derivation is only valid for qua- 
sistatic response above z c . 

To explain data such as in Fig. 1, we must identify a 
length scale £±(Az,w) valid not only above z c , but also 
below and at finite frequency w (or rate s). Here we de- 
termine £± and show that it is the size of a finite system in 
which elastic storage and viscous loss balance. We also 
provide scaling arguments showing that response func- 
tions such as J(s) and the complex shear modulus G*(lo) 
are controlled by £±. Moreover, their rate and frequency 
dependence matches those found in jammed solids [5], 
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FIG. 1: Numerical measurement of creep compliance J(s) 
in damped random spring networks for varying distance to 
isostaticity Az = z — z c . For long time (low s), networks 
below isostaticity begin to flow, J ~ 1/s (dashed line), while 
networks above isostaticity achieve a finite strain. 



biopolymer networks [9], and athermal suspensions [6j. 

Random spring networks. — The critical coordination 
of central force networks with dimensionality d is z c = 2d 
[1] . We generate networks near z c according to the proto- 
col of Ref. [5] , which begins from an initial network with 
high coordination and randomly removes springs to reach 
a targeted z. The initial disordered network is derived 
from a numerically generated sphere packing. During di- 
lution, springs are removed only if the affected nodes re- 
tain at least d+ 1 springs; we present numerical results for 
d = 2, though the scalings we predict are independent of 
d. This rule guarantees that all springs generically bear 
load in hyperstatic networks, and that even hypostatic 
networks are above connectivity percolation. 

A network's nodes are connected by springs at their 
rest length with stiffness k. When nodes undergo dis- 
placements {Ui}, the force from node j on node i, 



Jij — rCU 



ij n ij ' 



is along the unit vector n^ pointing from 



i to j and proportional to their relative normal motion 
u a = {Uj ~ ^i) ' n - To introduce damping, we take the 
network to be immersed in an affinely flowing Newtonian 
fluid: under a pure shear j(t), the fluid element at posi- 
tion x has displaced by ^{t)x at time t. A node's velocity 

relative to the fluid, or non-affine velocity U™, is damped 



by a viscous force F t = —bUf a : with damping coefficient 
b. We set the stiffness k, the damping b, and the average 
spring length all to unity. Equations of motion and nu- 
merical details are given in the Supplementary Material. 

Gradients in the solid and fluid phases. — Our scaling 
arguments are built on two simple relationships between 
relative normal motions u" and non-affine motions f/ na . 
We first develop these relationships before turning to £± . 

In a solid, zero frequency shear stretches and com- 
presses springs. The normal motion uL in spring (ij)is 
a discretization of V|| • U, the local gradient of the dis- 
placement field U along n. We assume the gradient is 
dominated by U na . For scaling purposes, we introduce a 
quantity A s , as yet unspecified, that relates the typical 
value of u" and the typical non-affine motion U na : 



U na 
^7 



(1) 



Because u" and U na are related by a gradient, X s is a 
length scale characterizing the solid (non-floppy) phase. 

A fluid possesses floppy modes and can perform zero 
frequency deformations without stretching springs. As 
there can be non-affine motion while u" is everywhere 
zero, Eq. (fTl) must break down. Elastic forces in a fluid 
are instead induced by flow. For oscillatory driving at 
finite frequency u> (or relaxation at rate s), nodes have 
finite velocity and experience viscous forces {-F»}. Since 
inertia is negligible for small ui, these viscous forces must 
be balanced by elastic forces {fij} in the springs. 

The force balance equations V ■ fij = Fj are a discrete 

counterpart of the continuum relation div <x = F, where 
F is the viscous force density and the stress tensor o 
is linear in the elastic forces. Thus viscous forces are 
related to gradients in the elastic forces, and there is an 
undetermined length scale A/ relating the typical elastic 
force / and viscous force F in the fluid (floppy) phase, 



F 



A/ 



(2) 



or equivalently u" ~ cuXfU na . 

In the following, an important consideration will be 
the complex shear modulus G*(u) = l/J(s)\ s=iU> = 
G'(lj) + iG"(ui). Its real and imaginary parts G' and 
G", known as the storage and loss moduli, are energy 
densities associated with elastic storage and viscous dis- 
sipation, respectively. Their scaling can be written 
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(A/C/ na ) uj 2 with floppy modes ,„, 

(t/ na /A s ) 2 without floppy modes, ^ ' 



where we have made use of Eqs. 11]) and (pi), and 

G" ~ F ■ U ~ (C/ na ) 2 u; , (4) 
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FIG. 2: (a) Three approaches to the critical point: along 
the z-axis from above and below, and along the frequency 
axis, (b) The diverging length scale £± has three regimes, 
distinguished by comparing the frequency uj to the diverging 
time scale r* ~ 1/Az 2 (dashed curve), (c) "Cutting out" 
a subsystem of size I. Springs crossing the dashed line are 
removed, (d) Creating a finite system by "freezing." Nodes 
in the shaded region are not allowed to move. 



both with and without floppy modes. Note that because 
G* (uj) and J(s) are related, for scaling purposes the in- 
verse time scales u> and s are interchangeable. 

Generalized isostatic length scale. — We now introduce 
a scaling argument to identify a dynamic critical length 
scale £±(Az,ui). The subscript indicates the sign of Az. 
In the conventional scenario for dynamic critical scal- 
ing near a critical point [13j . a diverging length scale 
displays three qualitatively distinct regimes (Fig. 2a,b). 
Two regimes persist at zero frequency and can be probed 
by approaching z c from above and below at ui = 0. The 
third or critical regime is inherently dynamical and can 
be probed by approaching z c along the frequency axis. 

Our approach is a generalization of the "cutting argu- 
ment," introduced in Refs. [2] and [TS], which identifies 
a diverging length scale £* ~ 1/Az associated with prox- 
imity to isostaticity. The derivation of t is appropriate 
only for hyperstatic systems at zero frequency. Because 
we generalize the cutting argument, our length scale re- 
covers I* in the quasistatic limit. 

We begin in the regime where £* obtains by approach- 
ing z c from above at w = 0. On long enough wavelengths 
it is reasonable to approximate a material as a contin- 
uum; £* is the length scale beyond which this approxi- 
mation is valid. One way to identify £* is to consider the 
properties of a finite portion of the infinite system. A 
sample with linear size £ that has the same properties as 
the infinite system is larger than £* . Alternatively, if the 
finite sample is qualitatively different - for example, if it 
no longer resists quasistatic shear - then £ < £*. Thus 
we must ask when such a difference first appears. 



So motivated, we imagine "cutting out" a sample of an 
infinite random spring network with excess coordination 
Az. To do this, we overlay a box of size £ on the infinite 
system and cut those springs that cross its boundary; see 
Fig. 2c. Cutting springs is destabilizing: cut enough and 
floppy modes will appear. When this happens the finite 
system must respond qualitatively differently to driving, 
as it cannot resist zero frequency shear. Therefore £* 
is the size of the box for which the first floppy mode 
appears. This occurs when the number of cut springs, 
proportional to the surface area £ d ~ 1 , exceeds the number 
of excess contacts, which is proportional to Az£ d . Hence 
(t) d - x ~ Az{£*) d , or £+(Az,0) =£* ~ 1/Az. 

We now approach z c from below. Clearly repeating the 
cutting approach will not work in a system that already 
has floppy modes. As cutting is in essence the introduc- 
tion of a free boundary, we now, instead, introduce a rigid 
boundary by "freezing" the position of nodes exterior to 
a box of linear size £ (see Fig. 2d). Freezing introduces 
stability: springs constrain floppy motion, and while each 
spring in the bulk is shared by two nodes, those springs 
connected at one end to a frozen node are not shared. 
Hence there are more constraints per degree of freedom in 
the finite system, and by choosing £ small enough we can 
remove all of its floppy modes. The counting is entirely 
analogous to the surface area-to-volume ratio described 
above and gives £_(Az,0) ~ l/(— Az). 

Lastly, we approach z c along the frequency axis. Recall 
that for the quasistatic systems, £±(Az,0) is the length 
scale for which a rigid system becomes floppy, or vice 
versa, under the influence of a boundary. As deforma- 
tions of random spring networks at finite frequency both 
store and dissipate energy, to proceed in analogy to the 
quasistatic case we should ask for what length scale stor- 
age and dissipation are comparable. A length scale can 
be imposed either by cutting or freezing. 

Let us first create a finite isostatic system by freezing. 
The frozen system has no floppy modes, so low frequency 
storage will dominate loss. We assume that the system 
organizes its response so as to minimize the stored energy 
G' ~ (U na J A s ) 2 . This is achieved when A s is as long as 
possible, X s ~ £ [J5]. By Eqs. (|3j and Q the ratio of 
loss to storage is then G"/G' ~ 1/£ 2 uj, which is of order 
unity when £ — £_(0,oj) ~ 1/uj 1 / 2 . We identify this 
length scale with £_ because it was reached via freezing. 
Alternatively, we can create a finite system by cutting an 
isostatic network. As cutting introduces floppy modes, 
low frequency loss will dominate storage. We assume that 
the response now minimizes the dissipated energy, G" ~ 
(J7 na ) 2 w ~ f 2 /X 2 f u, which is achieved if A/ is, likewise, 
identified with £ [T5]. The ratio G"/G' ~ £ 2 uj is again of 
order unity when £ = £ + (0,oj) ~ 1/uj 1 ' 2 . Interestingly, 
displacement correlations in normal modes of undamped 
hyperstatic networks also decay as 1/cj 1 / 2 [IT] . 

In summary, we have used cutting and freezing argu- 
ments to infer a dynamic critical length scale £±(Az,u>). 



The length scale has two branches, a solid branch for 
Az > and a fluid branch for Az < 0, which together 
can be expressed succinctly as 



£±(Az,w) 



±1/Az UJT* < I 



1/w 



1/2 



>1. 



(5) 



The crossover is controlled by the diverging time scale 
t* ~ 1/Az 2 , which follows from balancing the scaling 
of £± in different regimes. Because both branches of 
£± scale as 1/oj 1 ' 2 near z c , we can speak meaningfully 
of a single dynamic critical regime; see Fig. 2b. £± is 
frequency-independent outside the critical regime. 

Predicting response. — We now consider infinite or pe- 
riodic systems. Given £±, we can construct a scaling 
argument for the response of nearly isostatic networks. 

The key step is to identify the gradient scales A s and 
A/ in Eqs. (fTl) and pi) with £ + and £_, respectively. 
This ansatz, which is consistent with our approach to 
finite systems, can be directly tested in numerics. If £± 
scales as in Eq. pi) , then data from each branch will col- 
lapse when plotting Az a (v)\ /U na ) and Az a (f/F) versus 
uj/Az b for a = 1 and b = 2. As shown in Fig. 3a, there 
is excellent data collapse for a = 1.1(1) and b — 2.1(1), 
in good agreement with predictions. This confirms both 
the identification of A/ and A s with £± and the scaling 
arguments that led to £±. 

Having identified X s and A/, Eqs. pi) and Q relate 
the complex shear modulus to the typical magnitude of 
non-affine motion U na , which remains to be determined. 
As non-affinity is a form of fluctuation, it is reasonable 
to expect its divergence on approach to the critical point 
|13j . U na ~ £±. The exponent v can be determined from 
the dynamic viscosity in hypostatic networks - its scaling 
Vo ~ V( — Az) follows from a straightforward calculation 
(see Supplementary Material). Physically, the viscos- 
ity diverges because low frequency motions of hypostatic 
networks strongly project on their floppy modes; t]q is 
controlled by the floppy mode number density, which is 
proportional to z c — z. Recalling that 770 = lim^^o G" /u> 
and comparing to Eq. Q, we conclude that v = 1/2; this 
prediction is confirmed in the inset to Fig. 3a. 

With an expression for £± and the value of v, we have 
now determined the scaling of G' and G" . These can be 
expressed in an elegant and compact form by introducing 
elastic and viscous coefficients G± ~ l/£± and rj ± ~ £ ± . 
For hyperstatic networks Eqs. pi) and Q become 



G" - G+ and G" 



n+uj. 



(6) 



These are the moduli of the simplest possible viscoelas- 
tic solid, the Kelvin- Voigt solid - elastic and viscous el- 
ements with coefficients G+ and 77+ connected in paral- 
lel. Similarly, for hypostatic systems we recover the low 
frequency moduli of elements G_ and rj_ connected in 
series, i.e. a Maxwell fluid, 



G' - {n-uf/G- and G" - ??_w . 



(7) 
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FIG. 3: Scaling collapse of data from simulated networks. 
Legend as in Fig. 1. (a) The diverging length scale £±, deter- 
mined from ratios of mean non-affine motion to mean relative 
normal motion and mean elastic force to mean viscous force 
as a function of driving frequency u). The dashed curve has 
slope a/b. Inset: l/([/ na ) 2 vanishes linearly at z c in the zero 
frequency limit, (b) Scaling collapse of the compliance data 
in Fig. 1. The dashed curve has slope a' /&'. Inset: Elastic 
and viscous coefficients in the zero frequency limit. 



Thus a network's complex shear modulus is always de- 
scribed by an elastic element G± and a viscous dashpot 
■q±. Above isostaticity they are wired as a solid, below 
isostaticity they are wired as a fluid. 

Testing predictions. — In addition to the direct test of 
Fig. 3a, we provide two further tests of our predictions. 
The first probes the zero frequency limit and confirms the 
scaling of G± , which is given by lim w _>.o G above z c and 
lim w ^ (G"') 2 /G' below, and r)± = lim^^o G'/u, which 
is valid on both sides of z c . Fig. 3b (inset) shows that, 
as predicted, both G± and l/rj± grow linearly with Az 
on both sides of the transition. 

Our second test revisits the compliance data of Fig. 1; 
by collapsing J(s), we verify the scaling of £± on finite 
time scales. In hyperstatic networks the compliance is 
dominated by storage, J ~ 1/G+ ~ £+(Az,s), while in 
hypostatic systems it is controlled by loss, J ~ l/rjs ~ 
l/[£_(Az, s) s]. Therefore plotting Az° J versus s/Az b 
for a' — 1 and b' — 2 should collapse data on both sides of 
the critical point to a two-branched master curve with a 



critical regime J ~ 1/s 1 / 2 . We find excellent collapse for 
a' = 1.1(1) and b' = 2.0(1) (Fig. 3b), confirming these 
predictions. The 1/s 1 / 2 divergence implies that strain 
grows anomalously with the square root of time before 
achieving constant strain or strain rate. This is reminis- 
cent of Andrade creep in metals, although the Andrade 
creep exponent is 1/3 rather than 1/2 [2D] . 

Conclusions. — We have predicted, tested, and con- 
firmed the form of a diverging dynamic length scale near 
isostaticity. When combined with simple scaling argu- 
ments, this length scale correctly predicts the viscoelas- 
tic response of nearly isostatic networks. At frequencies 
low compared to 1/r*, £± is controlled by the distance to 
isostaticity, networks behave as simple viscoelastic solids 
or fluids, and their elastic and viscous moduli are set by 
Az. At higher frequencies there is a critical regime in 
which £± is determined by the time scale with which the 
system is driven, storage and loss are comparable, and 
the distinction between solid and fluid is blurred. 

Though we have treated spring networks, our predic- 
tions, in particular the rate-dependence of the compliance 
and shear modulus, match prior results for G* in jammed 
solids |8J and biopolymer networks 9] and the relaxation 
modulus in athermal suspensions [pj. This suggests the 
broader applicability of the scaling arguments presented 
here, which we expect can be extended to incorporate ad- 
ditional physics, including bending stiffness [5][TD], finite 
strain amplitude [H[T2J[5T], and steady flow [2"2"] . 

This work was supported by the Dutch Organization 
for Scientific Research (NWO). 



[1] 

[2] 



[10 

[11 
[12 

[13 

[14 



S. Alexander, Phys. Rep 296, 65 (1998). 

M. Thorpe, D. Jacobs, M. Chubynsky, and J. Phillips, 

Journal of Non-Crystalline Solids 266-269, Part 2, 859 

(2000), ISSN 0022-3093. 

C. Heussinger and E. Frey, Phys. Rev. Lett. 96, 017802 

(2006). 

M. Wyart, H. Liang, A. Kabla, and L. Mahadevan, Phys. 

Rev. Lett. 101, 215501 (2008). 

W. G. Ellenbroek, Z. Zeravcic, W. van Saarloos, and 

M. van Hecke, EPL (Europhysics Letters) 87, 34004 

(2009). 

T. Hatano, Phys. Rev. E 79, 050301 (2009). 

M. van Hecke, J. Phys. Cond. Matt. 22, 033101 (2010). 

B. P. Tighe, Phys. Rev. Lett. 107, 158303 (2011). 

C. P. Broedersz, M. Depken, N. Y. Yao, M. R. Pollak, 

D. A. Weitz, and F. C. MacKintosh, Phys. Rev. Lett. 
105, 238101 (2010). 

C. P. Broedersz, X. Mao, T. C. Lubensky, and F. C. 

MacKintosh, Nat Phys 7, 983 (2011). 

M. Wyart, EPL (Europhysics Letters) 89, 64001 (2010). 

M. Sheinman, C. P. Broedersz, and F. C. MacKintosh, 

Phys. Rev. E 85, 021801 (2012). 

B. I. Halperin and P. C. Hohenberg, Phys. Rev. 177, 952 

(1969). 

A. V. Tkachenko and T. A. Witten, Phys. Rev. E 60, 



[15 
[16 
[17 
[18 
[19 

po; 

[21 
[22 

[23 
[24 



687 (1999). 

M. Wyart, S. R. Nagel, and T. A. Witten, Euro- 
phys. Lett. 72, 486 (2005). 

L. E. Silbert, A. J. Liu, and S. R. Nagel, Phys. Rev. Lett. 
95, 098301 (2005). 

W. G. Ellenbroek, E. Somfai, M. van Hecke, and W. van 
Saarloos, Phys. Rev. Lett. 97, 258001 (2006). 
M. Mailman and B. Chakraborty, JSTAT 2011, L07002 
(2011). 

We expect that increasing A s increases £/ na and that in- 
creasing Xf decreases /, both of which are consistent with 
resulting predictions that are verified numerically. 
M.-C. Miguel, A. Vespignani, M. Zaiser, and S. Zapperi, 
Phys. Rev. Lett. 89, 165501 (2002). 

B. P. Tighe and J. E. S. Socolar, Phys. Rev. E 77, 031303 
(2008). 

B. P. Tighe, E. Woldhuis, J. J. C. Remmers, W. van Saar- 
loos, and M. van Hecke, Phys. Rev. Lett. 105, 088303 
(2010). 

C. E. Maloney and A. Lemaftre, Phys. Rev. E 74, 016118 
(2006). 

B. P. Tighe and T. J. H. Vlugt, J. Stat. Mech. p. P01015 
(2010). 



Supplementary Material. — For low frequency mo- 
tions inertia can be neglected. The overdamped equa- 
tions of motion can then be written |5] 



kq + iojBQ = anr . 



(8) 



Here a is the shear stress, f2 is the volume, and f is a unit 
vector along the strain coordinate. Q is a vector com- 
prising the dN displacement components and the shear 
strain 7. K and B are stiffness and damping matrices, 
respectively, that quantify the elastic and viscous inter- 
actions. They are matrices of second derivatives 



K„ 



d 2 V 



and Br, 



d 2 R 



(9) 



dQmdQn dQ m dQ 

where V is the elastic potential energy 

v = \Y, k ^l 3 ?' ( 10 ) 



and R is the Rayleigh dissipation function 

1, 



R 



\Hw 



'" 2 + -bj 2 V. 



(ii) 



The numerical results reported in Figs. 3 and 4 are 
for periodic networks of N = 256 nodes. For each co- 
ordination number we average over approximately 100 
networks. 

The viscoelastic response can be characterized by 
studying the generalized eigenvalues and -vectors of 
{K,B} [8 . The eigenvectors {U n } are evanescent with 
relaxation rates given by the eigenvalues {s n < 0}. 
Floppy modes have zero relaxation rate. Expanding Q 
in the eigenvectors and solving for the viscosity % yields 



N 

no 



E c «+ E 

(n|« n =0) (ra|s„<0) 



(12) 



The positive prefactor c„ = |U„ • f | 2 /(U„,BU„) is the 
strength of the coupling between eigenmode n and the 
imposed shear; it is positive and to very good approxi- 
mation independent of s n [8, 23J. It is straightforward 
to show that the sum over floppy modes is the dominant 
contribution to the viscosity, so that 



no 



N 



■Az' 



(13) 



where 7Vf m = —AzN/2 is the number of floppy modes 
in the hypostatic system. Thus the number density of 
floppy modes directly sets the viscosity in hypostatic net- 
works, similar to pressure fluctuations in jammed solids 
[24] . The dynamic viscosity, or equivalently the scaling 
of non-affine fluctuations, can also be calculated in qua- 
sistatic hyperstatic networks, where loss is subdominant 
[HIHIIZ]. One again finds 770 ~ 1/Az, consistent with 
the insets to Figs. 3a and b. Unlike the hypostatic case, 
however, these calculations require prior knowledge of ei- 
ther the static shear modulus or the form of the density 
of states. 
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